Introduction

  • Jeffrey D. Blume, PhD
    • School of Data Science, University of Virginia
  • Megan H. Murray, PhD
    • Department of Biostatistics, Vanderbilt University

Resources:

SGPV Package

Functions Included:

  • sgpvalue()
    • This function computes the second-generation p-value (SGPV) and its associated delta gaps, as introduced in Blume et al. (2018).
  • sgpower()
    • Calculate power and type I error values from significance testing based on second-generation p-values as the inferential metric.
  • plotsgpv()
    • This function displays user supplied interval estimates (support intervals, confidence intervals, credible intervals, etc.) according to its associated second-generation p-value ranking.
  • plotman()
    • This function displays a modified Manhattan-style plot colored according to second-generation p-value status.
  • fdrisk()
    • This function computes the false discovery risk (sometimes called the “empirical bayes FDR”) for a second-generation p-value of 0, or the false confirmation risk for a second-generation p-value of 1.

1 Part 2

1.1 Part 2a

1.1.1 Statistical Properties

3 inference outcomes:

  • Data are consistent with null \(p_\delta=1\)
  • Data are consistent with alternative \(p_\delta=0\)
  • Data are inconclusive \(0<p_\delta<1\)

2 underlying truths:

  • Null is true
  • Alternative is true

CIs can miss (bummer)

cidemo = function(rep=100,
                  mu=0.8,
                  n=200,
                  alpha=0.05,
                  ylb="Probability",
                  ltype=1,
                  lcol=4,
                  ltype1=1,
                  lwd1=2,
                  lcol1=2,
                  plotset=F,
                  hilim=1,
                  lolim=0.5,
                  bcol="lemonchiffon2",
                  xleg=45,
                  yleg=0.65,
                  idz.lo=0.75,
                  idz.hi=0.85,
                  zone=TRUE){

  cints=matrix(0,3,rep)
  
  cints[1,]=seq(1,rep,1)
  sigma=sqrt(mu*(1-mu))
  
  for (i in 1:rep) {
  
  dat=rbinom(n,1,prob=mu)
  
  ## Standard normal quantile
  za=abs(qnorm(alpha/2))
  
  ## upper CI limit
  cints[2,i]=mean(dat)+za*sigma/sqrt(n) 
  
  ## Lower CI limit
  cints[3,i]=mean(dat)-za*sigma/sqrt(n)
  }
  
  if (plotset==T){
    lolim=mu-5*sigma/sqrt(n)
    hilim=mu+5*sigma/sqrt(n)
  }
  
  if (sigma!=1){
    plot(cints[1,],cints[2,],ylim=c(lolim,hilim),type="n",
    main=paste("95% Confidence Intervals for the Probability\n Sample Size of ",n," per Endpoint"),
    ylab=ylb,xlab="Replication")
    }
    
  if (sigma==1){
    plot(cints[1,],cints[2,],ylim=c(lolim,hilim),type="n",
    main=paste("95% Confidence Intervals on the Standardized Effect Size \n Sample Size of ",n," per Endpoint"),
    ylab="Standardized Effect Size (X/sigma units)",xlab="Endpoint Identification Number")
    }

  if (zone==TRUE) {         
    polygon(c(0,rep,rep,0),c(idz.hi,idz.hi,idz.lo,idz.lo),density=200,col=bcol)
  }
            
  abline(h=mu)
  points(cints[1,],cints[2,],cex=0.6,pch=16,col=lcol)
  points(cints[1,],cints[3,],cex=0.6,pch=16,col=lcol)
  segments(cints[1,], cints[2,], cints[1,], cints[3,],lty=ltype,col=lcol)
  
  ct=0
  ct2=0
  for (i in 1:rep) {
  
    if (cints[2,i]<mu | cints[3,i]>mu){
        segments(cints[1,i], cints[2,i], 
                 cints[1,i],cints[3,i],
                 lty=ltype1,col=lcol1,lwd=lwd1) 
        points(cints[1,i],cints[2,i],cex=0.6,pch=16,col=lcol1)
        points(cints[1,i],cints[3,i],cex=0.6,pch=16,col=lcol1)
        ct=ct+1}
        
    if (cints[2,i]<idz.lo | cints[3,i]>idz.hi){
      ct2=ct2+1
    }
  }
  
  if (rep>100){
    abline(h=idz.hi,col=bcol,lty=2,lwd=1.5)
    abline(h=idz.lo,col=bcol,lty=2,lwd=1.5)
  }
  
  ## prob of missing indifference zone
  ## P(CI upper limit<idz.lo) +P(CI lower limit>idz.hi) 
  miss.idz=pnorm(-za+idz.lo*sqrt(n)/sigma)+pnorm(-za-idz.hi*sqrt(n)/sigma)
  
  if (zone==TRUE){
    legend(xleg,yleg,bty="n",
    c(paste("Missed True Value:             ",ct," / ",rep,sep=""),
    paste("Missed Indifference Zone:   ",ct2," / ",rep,sep="")))
  }
    
  if (zone!=TRUE){
    legend(xleg,yleg,bty="n",
    c(paste("Missed True Value:             ",
            ct," / ",rep,sep="")))
  }
  
}

set.seed(999)

cidemo(n=75,zone=FALSE)

cidemo(n=75,zone=FALSE)

cidemo(n=400,zone=FALSE)

cidemo(n=400,zone=FALSE)

cidemo(n=75)

cidemo(n=75)

cidemo(n=400)

cidemo(n=400)

prob.alt.sg = function(n, null.delta, 
                    theta0=0, 
                    theta=0,  
                    V=1, 
                    alpha=0.05){
  results = rep(NA, length(n))
  
  results = pnorm(sqrt(n)*(theta0-null.delta)/sqrt(V)-sqrt(n)*theta/sqrt(V)-qnorm(1-alpha/2))+
    pnorm(-1*sqrt(n)*(theta0+null.delta)/sqrt(V)+sqrt(n)*theta/sqrt(V)-qnorm(1-alpha/2))
  
  results
}

Errors and Plots

prob.null.sg = function(n, 
                        null.delta, 
                     theta0=0, 
                     theta=0,
                     V=1, 
                     alpha=0.05){
  results = rep(NA, length(n))
  
  if(length(null.delta)>1){
    con.large = which(null.delta>(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
    con.small = which(null.delta<=(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
    
    results[con.small] = rep(0, length(con.small))
    results[con.large] =
      pnorm(sqrt(n[con.large])*theta0+sqrt(n[con.large])*null.delta[con.large]/sqrt(V)-sqrt(n[con.large])*theta/sqrt(V)-qnorm(1-alpha/2))-
      pnorm(sqrt(n[con.large])*theta0-sqrt(n[con.large])*null.delta[con.large]/sqrt(V)-sqrt(n[con.large])*theta/sqrt(V)+qnorm(1-alpha/2))
  }else if(null.delta>(qnorm(1-alpha/2)*sqrt(V)/sqrt(n))){ 
    results =
      pnorm(sqrt(n)*theta0+sqrt(n)*null.delta/sqrt(V)-sqrt(n)*theta/sqrt(V)-qnorm(1-alpha/2))-
      pnorm(sqrt(n)*theta0-sqrt(n)*null.delta/sqrt(V)-sqrt(n)*theta/sqrt(V)+qnorm(1-alpha/2))
  }else{
    results = rep(0, length(theta))
  }
  results
}
prob.incon.sg = function(n, null.delta, 
                    theta0=0, 
                    theta=0,  
                    V=1,
                    alpha=0.05){
  results = rep(NA, length(n))
  
  con.large = which(null.delta>(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
  con.small = which(null.delta<=(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
  
  results = 1-(prob.alt.sg(n, null.delta,
                        theta0, 
                        theta, 
                        V,
                        alpha)+
                 prob.null.sg(n, null.delta,
                           theta0, 
                           theta,
                           V, 
                           alpha))
  
  results
}
theta.1 = seq(-3,3,by=0.01)

plot(theta.1, prob.alt.sg(n=50,
            null.delta = 0,
            theta=theta.1), type="l",
     ylim=c(0,1),
     xlab="Alternative",
     ylab="Probability",
     main="Prob of being consistent with alternative")
lines(theta.1, prob.alt.sg(n=50,
            null.delta = 0.5,
            theta=theta.1),
      col="green")
lines(theta.1, prob.alt.sg(n=50,
            null.delta = 1,
            theta=theta.1),
      col="blue")
abline(h=0.05, lty=2)
legend("bottomright",
       legend=c("Delta=0", "Delta=0.5", "Delta=1"),
       col=c("black", "green","blue"), lty=1)

theta.1 = seq(-3,3,by=0.01)

plot(theta.1, prob.null.sg(n=50,
            null.delta = 0,
            theta=theta.1), type="l",
     ylim=c(0,1),
     xlab="Alternative",
     ylab="Probability",
     main="Prob of being consistent with null")
lines(theta.1, prob.null.sg(n=50,
            null.delta = 0.5,
            theta=theta.1),
      col="green")
lines(theta.1, prob.null.sg(n=50,
            null.delta = 1,
            theta=theta.1),
      col="blue")
abline(h=0.05, lty=2)
legend("bottomright",
       legend=c("Delta=0", "Delta=0.5", "Delta=1"),
       col=c("black", "green","blue"), lty=1)

theta.1 = seq(-3,3,by=0.01)

plot(theta.1, prob.incon.sg(n=50,
            null.delta = 0,
            theta=theta.1), type="l",
     ylim=c(0,1),
     xlab="Alternative",
     ylab="Probability",
     main="Prob of being inconclusive")
lines(theta.1, prob.incon.sg(n=50,
            null.delta = 0.5,
            theta=theta.1),
      col="green")
lines(theta.1, prob.incon.sg(n=50,
            null.delta = 1,
            theta=theta.1),
      col="blue")
abline(h=0.05, lty=2)
legend("bottomright",
       legend=c("Delta=0", "Delta=0.5", "Delta=1"),
       col=c("black", "green","blue"), lty=1)

1.2 Part 2b

1.2.1 Equivalence tests: TOST

  • Establish bioequivalence between data and an established range
  • Example: A pharmaceutical company tests for drug approval by comparing new drug’s performance to an approved drug’s performance
  • Tests are ordinary, one-sided, α-level t-tests
  • Flips the null and alternative hypotheses
n=6
dat = rnorm(n, 0, 1)
a = t.test(dat)$conf.int[1]
b =  t.test(dat)$conf.int[2]
  
x_bar = mean(dat)
S = sd(dat)
Z = (b-x_bar)/(S/sqrt(n))
  
theta_p = 0.5*0.75
theta_m = -0.5*0.75
  
out_TOST = tsum_TOST(m1=x_bar, 
            mu=0,       
            sd1=S, 
            n1=n, 
            low_eqbound=theta_m, 
            high_eqbound=theta_p, 
            alpha=0.05,
            eqbound_type = "raw")

out_TOST$TOST
##                     t        SE df   p.value
## t-test      0.1645490 0.4629848  5 0.8757443
## TOST Lower  0.9745108 0.4629848  5 0.1872798
## TOST Upper -0.6454127 0.4629848  5 0.2735524
max(out_TOST$TOST$p.value[2], out_TOST$TOST$p.value[3])
## [1] 0.2735524
out_sgpv = sgpvalue(est.lo=a, 
                       est.hi=b, 
                       null.lo=theta_m, 
                       null.hi=theta_p)

out_sgpv 
## $p.delta
## [1] 0.5
## 
## $delta.gap
## [1] NA
set.seed(999)
iter=500
n=6
results = as.data.frame(matrix(ncol=3, nrow=iter))
colnames(results) = c("sgpv", "tost", "Case")
half.case = NULL
half.case2 = NULL

for(i in 1:iter){
  dat = rnorm(n, 0, 1)
  a = t.test(dat)$conf.int[1]
  b =  t.test(dat)$conf.int[2]
  x_bar = mean(dat)
  S = sd(dat)
  Z = (b-x_bar)/(S/sqrt(n))
  
  theta_p = 0.5*0.75
  theta_m = -0.5*0.75
  
  out_TOST = tsum_TOST(m1=x_bar, 
            mu=0,       
            sd1=S, 
            n1=n, 
            low_eqbound=theta_m, 
            high_eqbound=theta_p, 
            alpha=0.05,
            eqbound_type = "raw")
  
  out_sgpv = sgpvalue(est.lo=a, 
                       est.hi=b, 
                       null.lo=theta_m, 
                       null.hi=theta_p)$p.delta
  
  results[i,1] = out_sgpv 
  results[i,2] = pmax(out_TOST$TOST$p.value[2],out_TOST$TOST$p.value[3]) 
  
  if(theta_m<a & theta_p>a & theta_p<b){
    # Case 4 some overlap
    results[i,3] = "\n Case 4: \nPartial Overlap"
  }else if(theta_m<b & theta_p>b & theta_m>a){
    # Case 3 some overlap
    results[i,3] = "\n Case 3: \nPartial Overlap"
  }else if(theta_m>b | theta_p<a){
    # Case 1 no overlap
    results[i,3] = "\n Case 1: \nNo Overlap"
  }else if((theta_m>a & theta_p<b)|(theta_m<=a & theta_p>=b)){
    # Case 2 complete overlap
    results[i,3] = "\n Case 2: \nComplete Overlap"
  }
}


library(ggplot2)
scaleFUN <- function(x) sprintf("%.1f", x)

g = ggplot(results, aes(sgpv, tost))+
  xlab("SGPVs") +
  ylab("Equivalence p-values") +
  geom_point(aes(alpha=0.2)) +
  geom_abline(intercept = 1, slope = -1,linetype='dashed', alpha=0.6)+
  xlim(c(0,1))+
  scale_x_continuous(breaks= c( 0,0.120001,0.880001,1),labels=c("0", "~0.1", "~0.9", "1"), minor_breaks = c(0.12,0.88),
                      sec.axis=sec_axis(~./1, name="",breaks= c( 0.04,0.51,0.98), labels=c("Consistent \n with Alternative", "\n Inconclusive", "Consistent \n with Null")))+ 
  scale_y_continuous(breaks= c(0,0.07,1),  labels=c("0", expression(alpha), "1"),minor_breaks = c(0.070001),
                     sec.axis=sec_axis(~./1, name="", breaks= c(0.005,0.61), labels=c("Consistent \n with Null","Inconclusive")))+
  annotate(geom="text", x=0, y=0.9, label="D",
              size=5,
           family = "Courier-Bold")+
  annotate(geom="text", x=0.85, y=0.98, label="E",
              size=5,
           family = "Courier-Bold")+
  annotate(geom="text", x=0.91, y=0.98, label="F",
              size=5,
           family = "Courier-Bold")+
  annotate(geom="text", x=0, y=0, label="H",
              size=5,
           family = "Courier-Bold")+
  annotate(geom="text", x=0.15, y=0, label="I",
              size=5,
           family = "Courier-Bold")+
  annotate(geom="text", x=0.91, y=0, label="J",
              size=5,
           family = "Courier-Bold")+
  theme_bw()+
  theme(axis.text.y =  element_text(size=10, hjust=-0.5),
        axis.text.y.right =  element_text(size=10, hjust=0.7),
         axis.text.x.top =  element_text(size=10, vjust=2),
        axis.text.x =  element_text(size=10, vjust=-0.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_line(size = 0.8,colour="darkgrey"))+ 
    guides(alpha=FALSE)
g

ggplot(results, aes(sgpv, tost, col=Case, pch=Case))+
  geom_point(aes(), cex=1.5) +
  # geom_jitter(width = 0.02, height = 0.02)+
  # ggtitle(paste0("Derived Sample Size= ",n," Iterations= ", iter))+
  xlab("SGPV") +
  ylab("TOST Eq PV") +
  # facet_wrap(~ Case) +
  geom_abline(intercept = 1, slope = -1,linetype='dashed', alpha=0.6)+
  scale_color_manual(values=c("#0000B9",
                              "#00A300",
                              "#FD49A0",
                              "#FFA500",
                              "black"))+
  scale_shape_manual(values=c(16,21,17,24,21))+
    guides(alpha=FALSE)+
  scale_y_continuous(labels=scaleFUN)+ 
  scale_x_continuous(labels=scaleFUN)+
  guides(color = guide_legend(byrow = TRUE))+
   theme_bw()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.spacing.y = unit(0.3, 'cm'))

1.3 Part 2c

1.3.1 SGPV False Discovery Risk

  • False Discovery Rate (FDR)
    • SGPV = 0
  • False Confirmation Rate (FCR)
    • SGPV = 1
  • sgpv::fdrisk()
    • This function computes the false discovery risk (sometimes called the “empirical bayes FDR”) for a second-generation p-value of 0, or the false confirmation risk for a second-generation p-value of 1.
###### FDR rates
# false discovery risk with 95% confidence level
fdrisk(sgpval = 0,  
       null.lo = log(1/1.1), 
       null.hi = log(1.1),  
       std.err = 0.8,  
       null.weights = 'Uniform',  
       null.space = c(log(1/1.1), log(1.1)),  
       alt.weights = 'Uniform',  
       alt.space = 2 + c(-1,1)*qnorm(1-0.05/2)*0.8,  
       interval.type = 'confidence',  
       interval.level = 0.05)
## [1] 0.05949861
## with truncated normal weighting distribution
fdrisk(sgpval = 0,  
       null.lo = log(1/1.1), 
       null.hi = log(1.1),  
       std.err = 0.8,  
       null.weights = 'Point', 
       null.space = 0,  
       alt.weights = 'TruncNormal',  
       alt.space = 2 + c(-1,1)*qnorm(1-0.041/2)*0.8,  
       interval.type = 'likelihood',  interval.level = 1/8)
## [1] 0.04902575
# false discovery risk with LSI and wider null hypothesis
fdrisk(sgpval = 0,  
       null.lo = log(1/1.5), 
       null.hi = log(1.5), 
       std.err = 0.8,  
       null.weights = 'Point', 
       null.space = 0,  
       alt.weights = 'Uniform', 
       alt.space = 2.5 + c(-1,1)*qnorm(1-0.041/2)*0.8,  
       interval.type = 'likelihood',  
       interval.level = 1/8)
## [1] 0.01688343
# false confirmation risk example
fdrisk(sgpval = 1,  
       null.lo = log(1/1.5), 
       null.hi = log(1.5), 
       std.err = 0.15,  
       null.weights = 'Uniform',  
       null.space = 0.01 + c(-1,1)*qnorm(1-0.041/2)*0.15, 
       alt.weights = 'Uniform', 
       alt.space = c(log(1.5), 1.25*log(1.5)), 
       interval.type = 'likelihood',  
       interval.level = 1/8)
## [1] 0.03059522
###
##
#

False Discovery and Confirmation Rates

power.pv=function(n=20,sd=1,alpha=0.05,mu.l=-3,mu.u=3,delta=0.3,mu.0=0,r=1){

  ## Computes power when delta is zero and non-zero
  ## Inputs:
  ## n is sample size 
  ## indifference zone is: mu.0 +/- delta
  ## alternative space is: mu in [mu.l, mu.u]
  ## data ~ N(mu,sd)
  ## r = P(H1)/P(H0) for simple H1 and H0
  
  mu=seq(mu.l,mu.u,0.01)
  
  ## Power (null is mu.0 ; alt is mu)
  pow1.d=pnorm(-sqrt(n)*(delta-(mu.0-mu))/sd-qnorm(1-alpha/2))
  pow2.d=pnorm(-sqrt(n)*(delta+(mu.0-mu))/sd-qnorm(1-alpha/2))
  power.d=pow1.d+pow2.d
  
  ## Power without indifference zone (delta=0)
  pow1.0=pnorm(-sqrt(n)*(0-(mu.0-mu))/sd-qnorm(1-alpha/2))
  pow2.0=pnorm(-sqrt(n)*(0+(mu.0-mu))/sd-qnorm(1-alpha/2))
  power.0=pow1.0+pow2.0
  
  ## Type I error (null is mu.0 ; alt is mu.0)
  alpha.d=2*pnorm(-sqrt(n)*(delta)/sd-qnorm(1-alpha/2))
  
  ## Type I error without indifference zone (delta=0)
  alpha.0=2*pnorm(-sqrt(n)*(0)/sd-qnorm(1-alpha/2))
  
  ## P(sgpv=1|H0) not quite alpha b/c rules out inconclusive
  gam1.0=pnorm(sqrt(n)*(delta)/sd-qnorm(1-alpha/2))
  gam2.0=pnorm(sqrt(n)*(-delta)/sd+qnorm(1-alpha/2))
  gam.0=(gam1.0-gam2.0)*(delta>qnorm(1-alpha/2)*sd/sqrt(n))
  
  ##
  gam1.a=pnorm(sqrt(n)*(delta+mu.0-mu)/sd-qnorm(1-alpha/2))
  gam2.a=pnorm(sqrt(n)*(mu.0-delta-mu)/sd+qnorm(1-alpha/2))
  gam.a=(gam1.a-gam2.a)*(delta>qnorm(1-alpha/2)*sd/sqrt(n))
  
  ## False discovery rates (fdr1=FDR ; fdr0=FNDR)
  fdr.d = 1 / (1 + power.d/alpha.d * r)
  fndr.0 = 1 / (1 + (1-alpha.0)/(1-power.0) * (1/r))
  
  fdr.0 = 1 / (1 + power.0/alpha.0 * r)
  fndr.d = 1 / (1 + gam.0/gam.a * (1/r))
  
  ret=cbind(mu,alpha.d,power.d,fdr.d,fndr.d,
                alpha.0,power.0,fdr.0,fndr.0,gam.0,gam.a)
  colnames(ret)=c('mu','alpha.d','power.d','fdr.d','fndr.d',
                        'alpha.0','power.0','fdr.0','fndr.0',
                            'gam.0','gam.a')
  ret
}
par(mfrow=c(1,1))
prev.ratio=1
set.delta=0.5
set.n=16

all=power.pv(delta=set.delta,n=set.n,r=prev.ratio,alpha=0.05)

plot(all[,'mu'],all[,'fdr.0'],type="n",ylab="Probability",xlab="Alternative in SD units",
     ylim=c(0,prev.ratio/(prev.ratio+1)),xlim=c(-3,3))
title(main="False discovery and confirmation rates")

lines(all[,'mu'],all[,'fdr.0'],col="firebrick",lwd=1,lty=2)
lines(all[,'mu'],all[,'fndr.0'],col="dodgerblue",lwd=1,lty=2)

lines(all[,'mu'],all[,'fdr.d'],col="firebrick",lwd=2)
lines(all[,'mu'],all[,'fndr.d'],col="dodgerblue",lwd=2)

legend("topleft",c("",
                   expression(paste("P-value < ",0.05,sep="")),
                   expression(paste(P[delta]," = 0",sep="")),
                   "", "",
                   expression(paste("P-value > ",0.05,sep="")),
                   expression(paste(P[delta]," = 1",sep=""))),
       col=c("","firebrick","firebrick","","",
             "dodgerblue","dodgerblue"),
       lty=c(NA,2,1,NA,NA,2,1),lwd=c(NA,1,2,NA,NA,1,2),bty="n")

text(-2.875,0.5,"FDR")
text(-2.875,0.42,"FCR")

axis(side = 1, at = c(-set.delta,set.delta),cex.axis=0.8,col.ticks="black",lwd.ticks=1.5,tcl=0.6,labels=FALSE)
axis(side = 3, at = c(-set.delta,set.delta),cex.axis=0.8,col.ticks="black",lwd.ticks=1.5,tcl=0.6,labels=FALSE)

text(0,-0.0125,expression(paste(H[0]," Zone",sep="")),cex=0.8)
text(0,prev.ratio/(prev.ratio+1)*1.015,expression(paste(H[0]," Zone",sep="")),cex=0.8)

text(2.5,0.05/(prev.ratio+0.05)+.01,"alpha/(alpha+r)")

1.3.2 FDRestimation

Download the package from:

p.fdr.output = p.fdr(leukstats$p.value)

head(p.fdr.output$`Results Matrix`)
##     BH FDRs Adjusted p-values Raw p-values
## 1 0.1892484         0.1892484   0.05615324
## 2 0.9762477         0.9759186   0.93296848
## 3 0.6970430         0.6969894   0.47936375
## 4 0.2141365         0.2141365   0.06684255
## 5 0.9160462         0.9156422   0.81439179
## 6 0.6859496         0.6859383   0.46673059
summary(p.fdr.output)
## 
## Call:
## p.fdr(pvalues = leukstats$p.value)
## 
## Number of tests: 7128
## Raw p-value Range: [0, 0.99964]
## 
## Adjustment Method: BH
## False Discovery Rate Range: [0, 1]
## Findings at 0.05 level:
##    Significant (Reject): 1233
##    Inconclusive (Fail to Reject): 5895
## ---
## Estimated/Assumed Null Proportion (pi0): 1

Plots

plot(p.fdr.output)

plot(p.fdr.output, xlim=c(1100,1600), ylim=c(0, 0.2))

which(-1*(p.fdr.output$`Results Matrix`$`Adjusted p-values`)+p.fdr.output$fdrs>0.001)
## integer(0)
#Benjamini-Yeukatelli
p.fdr.output = p.fdr(leukstats$p.value, adjust.method = "BY")

head(p.fdr.output$`Results Matrix`)
##     BY FDRs BY Adjusted p-values Raw p-values   BH FDRs
## 1 0.1892484            0.1892484   0.05615324 0.1892484
## 2 1.0000000            1.0000000   0.93296848 0.9762477
## 3 1.0000000            1.0000000   0.47936375 0.6970430
## 4 0.4461177            0.4461177   0.06684255 0.2141365
## 5 1.0000000            1.0000000   0.81439179 0.9160462
## 6 1.0000000            1.0000000   0.46673059 0.6859496
plot(p.fdr.output)

plot(p.fdr.output, xlim=c(2000,2400), legend.on = FALSE)

Simple Example

pvalues=c(0.005,0.049,0.05,0.051,0.7)
zvalues=qnorm(pvalues/2, lower.tail = FALSE)

p.fdr.output = p.fdr(pvalues)

adj.pvalues= p.fdr.output$`Results Matrix`$`Adjusted p-values`
adj.fdrs= p.fdr.output$fdrs

single.fdr = c(p.fdr(pvalues=pvalues[1],zvalue=zvalues[1]),
               p.fdr(pvalues=pvalues[2],zvalue=zvalues[2]),
               p.fdr(pvalues=pvalues[3],zvalue=zvalues[3]),
               p.fdr(pvalues=pvalues[4],zvalue=zvalues[4]),
               p.fdr(pvalues=pvalues[5],zvalue=zvalues[5]))
df = data.frame("Raw p-values"= pvalues, 
                "Z-values" = zvalues,
                "BH Adj p-values" = adj.pvalues,
                "BH FDRs" = adj.fdrs,
                "Single lower bound FDRs" = single.fdr)

colnames(df) = c("Raw p-values", 
                "Z-values",
                "BH Adj p-values"  ,
                "BH FDRs" ,
                "Single lower bound FDRs" )
                
kable(round(df,3))%>%
  kable_styling(c("striped", "bordered"))
Raw p-values Z-values BH Adj p-values BH FDRs Single lower bound FDRs
0.005 2.807 0.025 0.025 0.019
0.049 1.969 0.064 0.122 0.126
0.050 1.960 0.064 0.083 0.128
0.051 1.951 0.064 0.064 0.130
0.700 0.385 0.700 0.700 0.481

Compare to p.adjust

set.seed(999)
pi0 <- 0.8
pi1 <- 1-pi0
n <- 10
n.0 <- ceiling(n*pi0)
n.1 <- n-n.0

sim.data <- c(rnorm(n.1,2,1),rnorm(n.0,0,1))
sim.data.p <- 2*pnorm(-abs(sim.data))

p.adjust.output = p.adjust(sim.data.p, method="fdr")

fdr.output = p.fdr(pvalues=sim.data.p, adjust.method="BH")

head(data.frame("Raw p-values"= sim.data.p,"p.adjust FDRs"=p.adjust.output,"p.fdr adj p-values"=fdr.output$`Results Matrix`$`Adjusted p-values`,"p.fdr FDRs"=fdr.output$fdrs))
##   Raw.p.values p.adjust.FDRs p.fdr.adj.p.values p.fdr.FDRs
## 1   0.08574923     0.4287462          0.4287462  0.4287462
## 2   0.49180527     0.7025790          0.7025790  0.7025790
## 3   0.42650649     0.7025790          0.7025790  0.7108441
## 4   0.78710602     0.7871060          0.7871060  0.7871060
## 5   0.78154483     0.7871060          0.7871060  0.8683831
## 6   0.57137764     0.7142221          0.7142221  0.7142221
plot(rank(sim.data.p, ties="random"), 
     fdr.output$`Results Matrix`$`Adjusted p-values`, 
     pch=17, 
     col="dodgerblue", 
     cex=2, 
     ylim=c(0,1),
     main="Comparison Plot", 
     xlab="Rank of p-values",
     ylab="")
points(rank(sim.data.p, ties="random"),p.adjust.output, pch=20, col="pink")
points(rank(sim.data.p, ties="random"),fdr.output$fdrs, pch=20, col="firebrick")
legend("bottomright", c("p.adjust FDRs", "p.fdr FDRs","p.fdr Adjusted p-values"), col=c("pink", "firebrick","dodgerblue"), pch=c(20,20,17))

Different Null proportions

get.pi0(leukstats$p.value, estim.method = "last.hist")
## [1] 0.6313131
get.pi0(leukstats$p.value, estim.method = "storey")
## [1] 0.6020751
get.pi0(leukstats$p.value, estim.method = "set.pi0", set.pi0=0.8)
## [1] 0.8
set.seed(999)
pi0 <- 0.8
pi1 <- 1-pi0
n <- 100
n.0 <- ceiling(n*pi0)
n.1 <- n-n.0

sim.data <- c(rnorm(n.1,2,1),rnorm(n.0,0,1))
sim.data.p <- 2*pnorm(-abs(sim.data))

fdr.output = p.fdr(pvalues=sim.data.p, adjust.method="BH")

plot(fdr.output)

fdr.output = p.fdr(pvalues=sim.data.p, adjust.method="BH", set.pi0=0.8)

plot(fdr.output)

get.pi0(sim.data.p, estim.method = "set.pi0", set.pi0=0.7)
## [1] 0.7
get.pi0(sim.data.p, estim.method = "last.hist")
## [1] 1
get.pi0(sim.data.p, estim.method = "storey")
## [1] 1
pi0 <- 0.8
pi1 <- 1-pi0
n <- 10000
n.0 <- ceiling(n*pi0)
n.1 <- n-n.0

sim.data <- c(rnorm(n.1,2,1),rnorm(n.0,0,1))
sim.data.p <- 2*pnorm(-abs(sim.data))

get.pi0(sim.data.p, estim.method = "set.pi0", set.pi0=0.7)
## [1] 0.7
get.pi0(sim.data.p, estim.method = "last.hist")
## [1] 0.796
get.pi0(sim.data.p, estim.method = "storey")
## [1] 0.8051147
Internal Code: Last Histogram Height Method
try.hist = hist(sim.data.p, breaks="scott")

try.mids = try.hist$mids
try.count = try.hist$counts

if(tail(try.mids,1)<0.5){
  pi0.hat=0
}else{
  pi0.hat = min(tail(try.count,1)*length(try.mids)/sum(try.count),1)
}

pi0.hat
## [1] 0.796

1.4 Part 2d

1.4.1 Study Planning

Purpose: Evaluate how different techniques of setting the indifference zone effect the SGPV study properties.

What options does a collaborator have when they are uncertain of the indifference zone?

Fixed Intervals

n= seq(1,250, by=0.1)
V=1
alpha=0.05

null.delta0 = rep(0.5, length(n))
null.delta1 = rep(0, length(n))
null.delta2 = rep(1.25, length(n))
null.delta3 = rep(1, length(n))
null.delta4 = rep(0.75, length(n))
null.delta5 = rep(0.25, length(n))

ci.bound = qnorm(1-alpha/2)*sqrt(V)/sqrt(n)
plot(n, null.delta1, 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,3),
     ylab=TeX(r'(Half Interval Width $ (\delta)$)'),
     xlab = "Sample Size (n)",
     main="")
lines(n, null.delta2,
      col="blue",
      lwd=2)
lines(n, null.delta3,
      col="purple",
      lwd=2)
lines(n, null.delta4,
      col="darkgreen",
      lwd=2)
lines(n, null.delta5,
      col="orange",
      lwd=2)
lines(n, null.delta0,
      col="hotpink",
      lwd=2)
lines(n, ci.bound,
      col="black",
      lty=2,
      lwd=2)
# abline(h=0.05, col="black", lty=2)
legend(130, 
       3,
        col=c("black",
             "red",
             "blue",
             "purple",
             "darkgreen",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(2,
             1,
             1,
             1,
             1,
             1,
             1),
       cex=0.65,
       lwd=2,
       legend=c(TeX(r'(Uncertainty Interval $\left(\sqrt{1/n}\right)$)'),
                "Point Null: Constant at 0",
                "SGPV: Constant at 1.25", 
                "SGPV: Constant at 1",
                "SGPV: Constant at 0.75",
                "SGPV: Constant at 0.5",
                "SGPV: Constant at 0.25"))

par(mfrow=c(1,2))
par(mar=c(2,2,2,2))

plot(n, prob.alt.sg(n, null.delta1,
                  theta=0), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,0.055),
     ylab="", 
     main=TeX(r'(Type 1 Error $\beta_0$)'))
lines(n, prob.alt.sg(n, null.delta2,
                  theta=0),
      col="blue",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta3,
                  theta=0),
      col="purple",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
                  theta=0),
      col="darkgreen",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
                  theta=0),
      col="orange",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
                  theta=0),
      col="hotpink",
      lwd=2)
abline(h=0.05, col="black", lty=2)

plot(n, prob.alt.sg(n, null.delta1,
                  theta=1), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Power $\beta_1$)'))
lines(n, prob.alt.sg(n, null.delta2,
                  theta=1),
      col="blue",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta3,
                  theta=1),
      col="purple",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
                  theta=1),
      col="darkgreen",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
                  theta=1),
      col="orange",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
                  theta=1),
      col="hotpink",
      lwd=2)
legend(70, 0.5,
       col=c("red",
             "blue",
             "purple",
             "darkgreen",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1,
             1),
       cex=0.7,
       lwd=2,
       legend=c("Point Null: Constant at 0",
                "SGPV: Constant at 1.25", 
                "SGPV: Constant at 1",
                "SGPV: Constant at 0.75",
                "SGPV: Constant at 0.5",
                "SGPV: Constant at 0.25"))

par(mfrow=c(1,1))
par(mar=c(4,4,2,2))

plot(n, prob.null.sg(n, null.delta1,
                  theta=0)+prob.alt.sg(n, null.delta1,
                  theta=1), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,2),
     ylab=TeX(r'( $\beta_1$ + $\omega_0$)'), 
     main= TeX(r'(Sum of Power and Probability of Null ($\beta_1$ +$\omega_0$))'))
lines(n, prob.null.sg(n, null.delta2,
                  theta=0)+prob.alt.sg(n, null.delta2,
                  theta=1),
      col="blue",
      lwd=2)
# lines(n, prob.null.sg(n, null.delta3,
#                   theta=0),
#       col="purple",
#       lwd=2)
lines(n, prob.null.sg(n, null.delta4,
                  theta=0)+prob.alt.sg(n, null.delta4,
                  theta=1),
      col="darkgreen",
      lwd=2)
lines(n, prob.null.sg(n, null.delta5,
                  theta=0)+prob.alt.sg(n, null.delta5,
                  theta=1),
      col="orange",
      lwd=2)
lines(n, prob.null.sg(n, null.delta0,
                  theta=0)+prob.alt.sg(n, null.delta0,
                  theta=1),
      col="hotpink",
      lwd=2)
legend(172, 0.7,
       col=c("red",
             "blue",
             "purple",
             "darkgreen",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1,
             1,
             2,
             2,
             2,
             2,
             2,
             2),
       cex=0.7,
       lwd=2,
       legend=c(TeX(r'( Constant at 0)'),
                TeX(r'(Constant at 1.25)'), 
                 TeX(r'( Constant at 1)'), 
                 TeX(r'( Constant at 0.75)'), 
                 TeX(r'(Constant at 0.5)'), 
                 TeX(r'( Constant at 0.25)')))

par(mfrow=c(1,2))
par(mar=c(2,2,4,2))

plot(n, prob.null.sg(n, null.delta1,
                  theta=0), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main= TeX(r'(Probability of Null $\omega_0$)'))
lines(n, prob.null.sg(n, null.delta2,
                  theta=0),
      col="blue",
      lwd=2)
lines(n, prob.null.sg(n, null.delta3,
                  theta=0),
      col="purple",
      lwd=2)
lines(n, prob.null.sg(n, null.delta4,
                  theta=0),
      col="darkgreen",
      lwd=2)
lines(n, prob.null.sg(n, null.delta5,
                  theta=0),
      col="orange",
      lwd=2)
lines(n, prob.null.sg(n, null.delta0,
                  theta=0),
      col="hotpink",
      lwd=2)

plot(n, prob.null.sg(n, null.delta1,
                  theta=1), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Probability of Null $\omega_1$)'))
lines(n, prob.null.sg(n, null.delta2,
                  theta=1),
      col="blue",
      lwd=2)
lines(n, prob.null.sg(n, null.delta3,
                  theta=1),
      col="purple",
      lwd=2)
lines(n, prob.null.sg(n, null.delta4,
                  theta=1),
      col="darkgreen",
      lwd=2)
lines(n, prob.null.sg(n, null.delta5,
                  theta=1),
      col="orange",
      lwd=2)
lines(n, prob.null.sg(n, null.delta0,
                  theta=1),
      col="hotpink",
      lwd=2)
legend(70, 0.5,
       col=c("red",
             "blue",
             "purple",
             "darkgreen",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1,
             1),
       cex=0.7,
       lwd=2,
       legend=c("Point Null: Constant at 0",
                "SGPV: Constant at 1.25", 
                "SGPV: Constant at 1",
                "SGPV: Constant at 0.75",
                "SGPV: Constant at 0.5",
                "SGPV: Constant at 0.25"))

par(mfrow=c(1,2))
par(mar=c(2,2,4,2))

plot(n,prob.incon.sg(n, null.delta1,
                  theta=0), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Probability of Inconclusive $\gamma_0$)'))
lines(n,prob.incon.sg(n, null.delta2,
                  theta=0),
      col="blue",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta3,
                  theta=0),
      col="purple",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
                  theta=0),
      col="darkgreen",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
                  theta=0),
      col="orange",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
                  theta=0),
      col="hotpink",
      lwd=2)

plot(n,prob.incon.sg(n, null.delta1,
                  theta=1), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Probability of Inconclusive $\gamma_1$)'))
lines(n,prob.incon.sg(n, null.delta2,
                  theta=1),
      col="blue",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta3,
                  theta=1),
      col="purple",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
                  theta=1),
      col="darkgreen",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
                  theta=1),
      col="orange",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
                  theta=1),
      col="hotpink",
      lwd=2)
legend(90, 0.8,
       col=c("red",
             "blue",
             "purple",
             "darkgreen",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1,
             1),
       cex=0.7,
       lwd=2,
       legend=c("Point Null: Constant at 0",
                "SGPV: Constant at 1.25", 
                "SGPV: Constant at 1",
                "SGPV: Constant at 0.75",
                "SGPV: Constant at 0.5",
                "SGPV: Constant at 0.25"))

Shrinking from 0.9 to 0.45

Let’s see how these look on a plot:

n2= seq(-10,300, by=0.1)
  
plot(n, ci.bound, 
     type="n",
     col="black",
     lty=2,
     lwd=2,
     ylim=c(-2,2.2),
     xlim=c(0,225),
     ylab="",
     xlab = "Sample Size (n)",
     main="",
     yaxt = "n")
# lines(n, -1*ci.bound,
#       col="black",
#       lty=2,
#       lwd=2)
lines(n, null.delta2,
      col="blue",
      lwd=2)
lines(n, -1*null.delta2,
      col="blue",
      lwd=2)
lines(n, null.delta4,
      col="hotpink",
      lwd=2)
lines(n, -1*null.delta4,
      col="hotpink",
      lwd=2)
lines(n, null.delta5,
      col="orange",
      lwd=2)
lines(n, -1*null.delta5,
      col="orange",
      lwd=2)
abline(h=k, col="black", lty=1,lwd=2)
abline(h=-1*k, col="black", lty=1,lwd=2)
abline(h=0, col="red", lty=1,lwd=2)
abline(h=1, col="red", lty=2,lwd=2)
polygon(c(n2, rev(n2)), c(rep(-1*k, length(n2)), rev(rep(-1*(k/2), length(n2)))),
        border = NA,
        col = rgb(210/633,210/633,210/633, alpha = 0.2))
polygon(c(n2, rev(n2)), c(rep((k/2), length(n2)), rev(rep(k, length(n2)))),
        border = NA,
        col = rgb(210/633,210/633,210/633, alpha = 0.2))
axis(2, at = c(-2,-1,-0.9,-0.45,0,
               2,1,0.9,0.45),
     labels = c(-2,-1,-0.9,-0.45,TeX(r'($\theta_0=0$)'),
               2,TeX(r'($\theta_1=1$)'),0.9,0.45),
     las=1)
legend(135, 
       2.25,
       col=c("red",
             "red",
             "black",
             "blue",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             2,
             1,
             1,
             1,
             1),
       cex=0.9,
       lwd=2,
       legend=c(TeX(r'(Null Hypothesis $\theta_0=0$)'),
                TeX(r'(Alternative Hypothesis $\theta_1=1$)'),
                "Original SGPV: Constant at 0.9", 
                TeX(r'($delta_1(n)\propto n^{-1/4}$)'),
                TeX(r'($delta_2(n)\propto n^{-1/3}$)'),
                 TeX(r'($delta_3(n)\propto n^{-1/2}$)')))

1.4.2 Probability Plots

par(mfrow=c(1,2))
par(mar=c(2,2,2,2))

plot(n, prob.alt.sg(n, null.delta1,
                  theta=0,
                  alpha=alpha.set), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,0.055),
     ylab="", 
     main=TeX(r'(Type 1 Error $\beta_0$)'))
lines(n, prob.alt.sg(n, null.delta2,
                  theta=0,
                  alpha=alpha.set),
      col="blue",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
                  theta=0,
                  alpha=alpha.set),
      col="black",
      lwd=2)
# lines(n, prob.alt.sg(n, null.delta3,
#                   theta=0),
#       col="hotpink",
#       lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
                  theta=0,
                  alpha=alpha.set),
      col="hotpink",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
                  theta=0,
                  alpha=alpha.set),
      col="orange",
      lwd=2)
abline(h=0.05, col="black", lty=2)

plot(n, prob.alt.sg(n, null.delta1,
                  theta=1,
                  alpha=alpha.set), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Power $\beta_1$)'))
lines(n, prob.alt.sg(n, null.delta2,
                  theta=1,
                  alpha=alpha.set),
      col="blue",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
                  theta=1,
                  alpha=alpha.set),
      col="black",
      lwd=2)
# lines(n, prob.alt.sg(n, null.delta3,
#                   theta=1),
#       col="hotpink",
#       lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
                  theta=1,
                  alpha=alpha.set),
      col="hotpink",
      lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
                  theta=1,
                  alpha=alpha.set),
      col="orange",
      lwd=2)
legend(70, 0.8,
       col=c("red",
             "black",
             "blue",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1),
       cex=0.7,
       lwd=2,
       legend=c("Point Null: Constant at 0",
                "Original SGPV: Constant at 0.9",
                TeX(r'($n^{-1/4}$)'),
                TeX(r'($n^{-1/3}$)'),
                TeX(r'($n^{-1/2}$)')))

par(mfrow=c(1,2))
par(mar=c(2,2,2,2))

plot(n, prob.null.sg(n, null.delta1,
                  theta=0,
                  alpha=alpha.set), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main= TeX(r'(Probability of Null $\omega_0$)'))
lines(n, prob.null.sg(n, null.delta2,
                  theta=0,
                  alpha=alpha.set),
      col="blue",
      lwd=2)
lines(n, prob.null.sg(n, null.delta0,
                  theta=0,
                  alpha=alpha.set),
      col="black",
      lwd=2)
# lines(n, prob.null.sg(n, null.delta3,
#                   theta=0),
#       col="hotpink",
#       lwd=2)
lines(n,prob.null.sg(n, null.delta4,
                  theta=0,
                  alpha=alpha.set),
      col="hotpink",
      lwd=2)
lines(n, prob.null.sg(n, null.delta5,
                  theta=0,
                  alpha=alpha.set),
      col="orange",
      lwd=2)

plot(n, prob.null.sg(n, null.delta1,
                  theta=1,
                  alpha=alpha.set), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Probability of Null $\omega_1$)'))
lines(n, prob.null.sg(n, null.delta2,
                  theta=1,
                  alpha=alpha.set),
      col="blue",
      lwd=2)
lines(n, prob.null.sg(n, null.delta0,
                  theta=1,
                  alpha=alpha.set),
      col="black",
      lwd=2)
# lines(n, prob.null.sg(n, null.delta3,
#                   theta=1),
#       col="hotpink",
#       lwd=2)
lines(n, prob.null.sg(n, null.delta4,
                  theta=1,
                  alpha=alpha.set),
      col="hotpink",
      lwd=2)
lines(n, prob.null.sg(n, null.delta5,
                  theta=1,
                  alpha=alpha.set),
      col="orange",
      lwd=2)
legend(70, 0.75,
       col=c("red",
             "black",
             "blue",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1),
       cex=0.7,
       lwd=2,
       legend=c("Point Null: Constant at 0",
                "Original SGPV: Constant at 0.9",
                TeX(r'($n^{-1/4}$)'),
                TeX(r'($n^{-1/3}$)'),
                TeX(r'($n^{-1/2}$)')))

par(mfrow=c(1,2))
par(mar=c(2,2,2,2))

plot(n,prob.incon.sg(n, null.delta1,
                  theta=0,
                 alpha=alpha.set), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Probability of Inconclusive $\gamma_0$)'))
lines(n,prob.incon.sg(n, null.delta2,
                  theta=0,
                  alpha=alpha.set),
      col="blue",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
                  theta=0,
                  alpha=alpha.set),
      col="black",
      lwd=2)
# lines(n,prob.incon.sg(n, null.delta3,
#                   theta=0),
#       col="hotpink",
#       lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
                  theta=0,
                  alpha=alpha.set),
      col="hotpink",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
                  theta=0,
                  alpha=alpha.set),
      col="orange",
      lwd=2)

plot(n,prob.incon.sg(n, null.delta1,
                  theta=1,
                 alpha=alpha.set), 
     type="l",
     col="red",
     lwd=2,
     ylim=c(0,1),
     ylab="", 
     main=TeX(r'(Probability of Inconclusive $\gamma_1$)'))
lines(n,prob.incon.sg(n, null.delta2,
                  theta=1,
                  alpha=alpha.set),
      col="blue",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
                  theta=1,
                  alpha=alpha.set),
      col="black",
      lwd=2)
# lines(n,prob.incon.sg(n, null.delta3,
#                   theta=1),
#       col="hotpink",
#       lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
                  theta=1,
                  alpha=alpha.set),
      col="hotpink",
      lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
                  theta=1,
                  alpha=alpha.set),
      col="orange",
      lwd=2)
legend(70, 0.55,
       col=c("red",
             "black",
             "blue",
             "hotpink",
             "orange"),
       box.lty=0,
       lty=c(1,
             1,
             1,
             1,
             1),
       cex=0.7,
       lwd=2,
       legend=c("Point Null: Constant at 0",
                "Original SGPV: Constant at 0.9",
                TeX(r'($n^{-1/4}$)'),
                TeX(r'($n^{-1/3}$)'),
                TeX(r'($n^{-1/2}$)')))

1.4.3 COVID Clinical Trail Example

  • Randomized 1,591 patients to ivermectin treatment or placebo
  • Mean time spent unwell was estimated using a longitudinal ordinal regression model; range was 0 to 14 days
  • Patients reported each day their symptoms and severity, health care visits, and medications.
  • “The difference in the amount of time spent feeling unwell with COVID was estimated to be 0.49 days in favor of ivermectin with a 95% credible interval of (0.15, 0.82).”
#Setup
theta0 = 0
point = 0.49 
data.lo = 0.15
data.hi = 0.82 
n=1591

data.sd=0.17

library(sgpv)

1 day difference

sgpvalue(est.lo = data.lo, 
         est.hi = data.hi,
         null.lo = -1,
         null.hi = 1)
## $p.delta
## [1] 1
## 
## $delta.gap
## [1] NA

Shrinking

k=1.5
V=data.sd^2

shrink.int=(k-(k/2))*sqrt(V)/(n^(1/2))+(k/2)

sgpvalue(est.lo = data.lo, 
         est.hi = data.hi,
         null.lo = -1*shrink.int,
         null.hi = shrink.int)
## $p.delta
## [1] 0.9002933
## 
## $delta.gap
## [1] NA
k=0.75

shrink.int=(k-(k/2))*sqrt(V)/n^(1/2)+(k/2)

sgpvalue(est.lo = data.lo, 
         est.hi = data.hi,
         null.lo = -1*shrink.int,
         null.hi = shrink.int)
## $p.delta
## [1] 0.3382063
## 
## $delta.gap
## [1] NA